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ABSTRACT 

We estimate the bispectrum of the Very Small Array data from the compact and ex- 
tended configuration observations released in December 2002, and compare our results 
to those obtained from Gaussian simulations. There is a slight excess of large bispec- 
trum values for two individual fields, but this does not appear when the fields are 
combined. Given our expected level of residual point sources, we do not expect these 
to be the source of the discrepancy. Using the compact configuration data, we put 
an upper limit of 5400 on the value of /nl 7 the non-linear coupling parameter, at 95 
per cent confidence. We test our bispectrum estimator using non-Gaussian simulations 
with a known bispectrum, and recover the input values. 

Key words: cosmology: observations - methods: data analysis - cosmic microwave 
background 



1 INTRODUCTION 

Currently-favoured cosmological theories predict that the 
primordial fluctuations in the cosmic microwave back- 
ground (CMB) obey Gaussian statistics to a high degree. 
The majority of inflationary scenarios imply a level of 
non-Gaussianity that is unlikely to be detectable by any 
forthcoming experiment ijAcauaviva et al1l2003t iMaldacenal 
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2003), although non-linear gravitational evolution may re- 
sult in detectable non- Gaussianity from the initially almost - 
Gaussian fluctuations feartolo. Matarrese fc Riottol Eooil) . 
Any convincing evidence for a departure from primordial 
Gaussianity would therefore play a very significant role 
in constraining theories of inflation. However, the dom- 
inant contributions to non-Gaussianity in the CMB are 
expected to come from secondary effects such as grav- 
itational lensing, reionization, the Sunyaev-Zel'dovich ef- 
fect, and from the local Universe. These effects are of 
varying significance, depending on the scale. Coupling be- 
tween lensing and the Sunyaev-Zel'dovich effect tends to 
dominate over primordial non-Gaussianity on small scales 
jGoldberg fc Spergell Il99sh . The effects of dust and gas 
clouds are very dependent on the region of sky observed. 
The Very Small Array (VSA) fields were carefully chosen to 
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minimise contamination from Galactic emission and bright 
radio sources, and the signal from Galactic foregrounds is 
thought to be much less than the primordial CMB signal, 
with the amount of contamination smaller on smaller scales 
jTavlor et al]l2003h . The recent results fr om the Wilkin- 
son M icrowave Anisotropy Probe ( WMAP; iKomatsu et alJ 
2003) have considerably tightened the constraints on primor- 
dial non-Gaussianity, whilst detecting a non-Gaussian signal 
arising from residual point sources. The non-Gaussian sig- 
nal that was found in the bispect rum of the COBE data 
iFerreira. Maeueiio fc Go rski 1998) has not been replicated 
in the WMAP data and is now believed to b e a result of sys- 
tematic errors dMaeueiio fc Medeirosll2003l) . The VSA has a 
dedicated point source subtractor, in order to remove all 
the sources above a certain minimum flux level, so that the 
residual contributio n from unsubtracte d sources is less than 
the flux sensitivity (iTavlor et al.ll2003l) . 



Whilst we may therefore have little basis for expecting 
to detect primordial non-Gaussianity in the VSA data, it 
is important to test the data nevertheless, if only to val- 
idate the assumption of Gaussianity which is made dur- 
ing the estimation of the power spectrum and its errors 
jHobson fc Maisingerll2002l). Currently the pu blished VSA 
data extend past I = 1400 (Grain ge et al.l 2003), m compar- 
ison with the current WMAP data which do not go beyond 
t = 900, so we are probing the fluctuations at higher reso- 
lution. Testing the data may also help to ascertain whether 
point source subtraction has been performed satisfactorily, 
and to see if there is any contamination by foreground 
sources. 



The VSA data have already been tested for non- 
Gaussianity using a variety of statistics in the map plane, 
and in the visibility plane by adopting a non- Gaussian like- 
lihood function. The results are presented in ISavage et alJ 
i2004h . Here, we use the bispectrum, the three-point statis- 
tic in the visibility plane. 

In Section |3 we discuss the statistics of CMB temper- 
ature fluctuations, and how the level of non-Gaussianity 
can be measured by the bispectrum. We then give a brief 
overview of the VSA in Section particularly with refer- 
ence to the point source subtraction technique, and present 
an approximate expression for the measured visibility three- 
point function, taking into account the convolution with the 
primary beam. An exact expression is derived in Appendix 
O In Section O we describe our method for estimating 
the bispectrum of the VSA data, which we have tested by 
producing non-Gaussian simulations as described in Section 
14.21 We then discuss a method for comparing the results to 
Gaussian simulations in Section ^. 3l and present our results. 
In Sections 14. 4l and l4.5l we consider the effect of point sources 
on the bispectrum, and look at the feasibility of detecting 
them in this way. Finally, in Section 14.61 we investigate the 
constraints that the VSA data are able to place on primor- 
dial non-Gaussianity. In the appendices we discuss optimal 
cubic bispectrum estimators, and the effect of the primary 
beam on the observed bispectrum from interferometric data. 



2 STATISTICS OF CMB TEMPERATURE 
FLUCTUATIONS 

2.1 Power spectrum 

Assuming full sky coverage, the temperature fluctuations 
of the CMB can be decomposed into spherical harmonics 
(Yem), and hence expressed as 

— (n) = ^ ai m Y lm (n) . (1) 

l,m 

Here, T = To, the mean temperature of the CMB. 
Rotational invariance demands that (p.t yrny a1 mo } = 
Ce 1 5e 1 e 2 Sm 1 m 2 , where the brackets denote the ensemble av- 
erage value. The values of the Ce represent the power spec- 
trum of the CMB, which has been the major focus of at- 
tention in CMB studies. The measured values of the Ce will 
differ slightly from the ensemble-averaged Ce due to instru- 
ment noise and cosmic variance; if the ae m are each drawn 
independently from a Gaussian distribution, with mean 
and variance Ce, then the measured Ce will be drawn from 
a \ 2 distribution, with an intrinsic variance equal to 2e°+i 
Therefore, particularly at low I, there will always be an un- 
certainty as to the true value of the Ce, independent of the 
technical difficulties of experimental measurement. 

2.2 Higher-order statistics 

If the temperature fluctuations of the CMB are Gaus- 
sian, the power spectrum completely describes the statistics. 
However, if there is a departure from Gaussianity, higher- 
order statistics are needed for a full description. In princi- 
ple there is an infinite number of ways in which the CMB 
could be non-Gaussian, and therefore the optimal statistic 
to use depends on the type of non-Gaussianity present. If 
we are looking for a particular signal we can use a statis- 
tic which is tailored for optimal detection of that signal. 
This method has been used to detect non- Gaussianity aris- 
ing fr om point sources in the WMAP data l|Komatsu et alJ 
2003). However, if we are not seeking a specific signature, 
then we require something which is fairly general. The nat- 
ural follow-on from the power spectrum is the bispectrum, 
which is the three-point function given by 

£1*2*3 = \ a eim 1 a e 2 m.2 a e 3m3 ) ■ 

The bispectrum gives a scale-dependent measure of skew- 
ncss bantc >s et al]l2003h . The ensemble-average bispectrum 
will be zero in the case of Gaussian fluctuations. However, 
for a given realisation, it will be non-zero, even neglecting 
instrumental noise and resolution, owing to cosmic variance. 

The ensemble-averaged bispectrum is constrained by 
the assumption of rotational invariance in a similar way to 
the power spectrum. This means that all the information is 
contained in the dependence of the bispectrum on the values 
of £, and not in its dependence on m. It can be shown that 

33/ V mi m% 7713 / 

1 In general the ae m are complex, and the variance of the real 
and imaginary parts is Ci/2, but since they satisfy the relation 
a lm = ( — l) ma J_ m they have 21+1 degrees of freedom. 



© 2004 RAS, MNRAS 000.IT1IT71 



Estimating the bispectrum of the Very Small Array data 3 



where the matrix re presents the Wigner 3-j symbol 
llRotenberg et al1ll959l) . Using the relation 
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we find that we can estimate the bispectrum as 
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If the CMB is Gaussian, then the cosmic variance of the 
bispectrum can be shown to be 

^Ws) (6) 
— C^CtzCtz (1 + 28e 1 £ 2 St 2 t 3 + S ei £ 2 + 5e 2 c a + 5e 3 c 1 ) 

The variance is twice as large if two Is are the same, or 
six times as large if all the Is are equal. This variance is 
enhanced if we only observe a portion of the sky (sample 
variance). In real experiments the measured signal will have 
a contribution arising from instrumental noise, which will 
be an additional source of variance. It is usually reasonable 
to assume that the noise is Gaussian, but it is conceivable 
that this could contribute to a non-zero bispectrum. 

In general, we can define any n-point function in har- 
monic space in a similar way to the power spectrum and 
bispe ctrum, and use these as a test for non-Gaussianity ijHiJ 
2002). However, definite detection of a signal arising from 
true non-Gaussianity at higher orders is progressively more 
unlikely as n increases. 



2.3 Flat-sky approximation 

If we are only observing a s mall patch of sky then we can use 
the fiat-sky approximation l)Hull200oT) . Instead of decompos- 
ing the temperature fluctuation into spherical harmonics, we 
instead use Fourier modes, so that the temperature fluctua- 
tions can now be expressed as 

AT f 

—p^- (x) — / a(u) exp(27riu ■ x)d u. (7) 

Translational and rotational invariance demand that 
(a{u)a{u'Y) = C(\u\)5' 2 (u - u'), 

and there is a straightforward correspondence with the full 
sky power spectrum at large £ : C(|u|) as Ce \i=2-k\u\- Simi- 
larly, assuming rotational and parity invariance, 

(a(u 1 )o(« 2 )o(« 3 )) = B(e 1 ,£ 2 ,(-3)S 2 (u 1 +u 2 + u 3 ), (8) 

with the correspondence li = 27r|tii|. Isotropy demands 
that the bispectrum is zero unless the three M-vectors sum 
to zero, so the bispectrum is specified by just two vec- 
tors, hence its name. Isotropy and parity invariance demand 
that the permutation-symmetric B{£\,£2,£3) depends only 
on the lengths of u\,U2 and 113. At large I we can write 
B(£i, £2, £3) ~ btuii-,, where ht, t 2 t a is the reduced bispec- 



trum iKomatsu fc SpergeJl2001T) 



Br 



(2ii + l) (24+1) (24+1) 
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The size of the VSA primary beam in the compact con- 
figuration implies that the flat-sky approximation is not en- 
tirely accurate. However, the correction factor is small, as 



described in lHobson fc Maisingerl ((2002) . The correction can 
be made be redefining the primary beam, and this is done 
in the power spectrum analysis. Our current bispectrum cal- 
culation does not directly take this effect into account, but 
we anticipate the associated error in the approximation to 
be small. 



3 THE VERY SMALL ARRAY 

The VSA is a 14-element interferometer situated on Mount 
Teide in Tenerife. It has the ability to make observations 
anywhere in the frequency range 26-36 GHz with an ob- 
serving bandwidth of 1.5 GHz. During the first observing 
season (September 2000 - September 2001) the antennae, 
of FWHM 4? 6, were arranged in a compact configuration, 
with a maximum baseline of 1.5 m, and observations were 
made at 34 GHz. In this configuration the VSA was sensitive 
to angular scales corresponding to £ ~ 150-900. The VSA 
was then reconfigured to form the extended array, with an- 
tennae of FWHM 2? 05 separated by a maximum distance of 
2.5 m, sensitive to the range £ ~ 300-1400. It has been mak- 
ing observations in this configuration since October 2001, at 
33-34 GHz. 

The compact array spent most of its time observing 
three separate regions of sky (labelled VSA1, VSA2 and 
VSA3). Within each region it made overlapping (mosaiced) 
observations of two or three fields (denoted by either no suf- 
fix, or t he suffixes A ; B or - OFF). Details of the fields are 
given in lTavlor et al.l l)2003|) . and the power sp ectrum cal- 
culate d from these observations is presented in lScott et alJ 
(2003). 

The extended array made mosaiced observations of 
smaller fields within the same three regions of the sky (la- 
belled with the suffixes E , F an d G), details of which can 
be found in lOrainge et all 120031. Here, we have studied the 
compact and extended array data that were used to com- 
pute t he CMB power spectrum presented in lGrainge et all 
(2003). Since completing the analysis here there has been a 
new relea se of results from further extended array observa- 
tions dPickinson et all2004h . Further details of the VSA ob - 
servational technique can be found in lWatson et alJ {2003). 

3.1 Interferometer measurements 

An interferometer samples the Fourier transform of the sky 
intensity multiplied by the primary beam. The fact that an 
interferometer measures directly in Fourier space makes in- 
terferometric data particularly suited to analysis in Fourier 
space. The visibility measured by the interferometer can be 
expressed as 



V{u) = / AI(x) A{x) e 27 ""'* d 2 x + N(u), 



(9) 



where AI(x) is the intensity fluctuation, A(x) the primary 
beam and N{u) the noise on baseline u. To a good approx- 
imation the beamshape can be modelled as Gaussi an with 
A(x) = exp(-|a;| 2 /2(T 2 ) ijHobson fc Maisingerll2002h . which 
is the shape that we have assumed when simulating the ob- 
servations. In Appendix[B]we describe how this leads to the 
approximate expression for the variance of the visibilities 



(V(u)V*(u))^na 2 f 2 C(u)+a 2 u , 



(10) 
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where 



4 CALCULATION OF THE BISPECTRUM 



f = T 



dB{v,T) 



dT 



(11) 



T=T 



is the conversion to specific intensity and a\ is the variance 
of the noise. Similarly, for the three-point function we obtain 



(V(tti)V(u3)V(«s)) 



f 3 ^na 2 B(u 1 ,u 2 ,u 3 ) 



(12) 



for ui + U2 + ti3 = 0. This approximation has been used to 
convert the measured three-point function to a bispectrum 
in units of /xK 3 for the purposes of plotting. In Appendix 
10 we derive an exact expression for a general case with a 
Gaussian beam which it is necessary to use when trying to 
quantify a true level of non-Gaussianity and considering a 
particular non-Gaussian signal. 

The sample time per visibility measurement by the VSA 
is typically 64 seconds. Each field was observed for ~100- 
200 hours, resulting in a very large number of individual 
visibility measurements. To reduce the data to a manageable 
level, the measurements are binned into square cells on the 
u v plane, using a maximum-l ikelihood method as described 
in lHobson fc Maisingerl J2002T) . Adjacent cells are still highly 
correlated as the cell width is chosen to be less than the 
width of the aperture function (the Fourier transform of the 
primary beam) so as not to lose information. 



4.1 Method of estimating the bispectrum 

Our starting point for estimation of the bispectrum of the 
VSA data is equation ® . The delta function tells us that we 
should be interested only in sets of three vectors in the uv 
plane that form a closed triangle. The lengths of the sides of 
the triangle give us the values of 1%, £2 and £3. It follows that 
we can make an estimate of B(£i,£2,£ 3 ) by finding from the 
data all possible triangles with these sides. 

There are two additional factors to consider: 

(i) The finite beamwidth means that the measured visibil- 
ity on baseline u is a convolution with the aperture function 
(the Fourier transform of the primary beam), as shown by 
equation 1B5II . so the visibilities are 'smeared' over a region 
of diameter a few wavelengths on the uv plane. Therefore it 
is not desirable to form separate bispectrum estimates using 
nearby visibility points. It is preferable to form individual 
estimates from 'bands' of points on the uv plane, else ad- 
jacent estimates will be highly correlated. This also means 
that we can slightly relax the requirement that the three 
vectors must form a closed triangle, as long as the deviation 
from a closed triangle is significantly less than the size of 
the aperture function. 

(ii) The quantity that we are interested in is the ensemble- 
averaged value of B(£i,£ 2 ,£ 3 ), however, since we are observ- 
ing only a small region of our Universe, we want to average 
over as many triangles as is reasonable in order to reduce 
the variance of the estimate. 



3.2 Source subtraction 

The VSA regions were carefully chosen in order to have 
relatively low dust contamination (using the dust maps of 
ISchlegel. Finkbeiner fc Davisl Il998h . low Galactic free-free 
and synchrotr on emission (as pred icted by the 408-MHz all- 
sky survey of lHaslam et aljf l982l and to avoid bright ra- 
dio sources l|TavloretHnTl2003F) . Choosing to observe at the 
higher end of the VSA fr equency range redu ces the signal 
from Galactic foregrounds ijTavlor et all2003T) . The compact 
array fields were chosen so as to contain no sources brighter 
than 0.5 Jy. In order to eliminate contamination from point 
sources, which would otherwise contribute to excess power, 
they are observed and directly subtracted from the data. 
The source subtraction strategy used was first to survey the 
fields with the Ryle telescope in Cambridge at 15 GHz prior 
to observation with the VSA. Then, simultaneous with VSA 
observations, a single-baseline interferometer, with a dish 
separation of 9 m, situated alongside the VSA, is used to 
monitor the sources identified in the 15-GHz survey. For the 
compact array, it was calculated that all sources brighter 
than 80 mjy needed to be subtracted. Therefore the survey 
with the Ryle telescope sought to identify all sources above 
20 mjy at 15 GHz. The power of the CMB is more sensi- 
tive to sources at higher values of £, and so for the extended 
array the fields were chosen, using the previ ous survey, to 
conta in no sources brighter than 100 mjy dGrainge et alJ 
2003). The sensitivity of the Ryle survey was increased to 
10 mjy, and all sources brighter than 20 mjy subtracted 
from the measured visibilities. 



In a general non-Gaussian scenario we would expect 
B(£i,£2, £3) to vary reasonably slowly with £; if it fluctuates 
too rapidly then the signal will be 'washed out' both by the 
convolution and by combining nearby points into one esti- 
mate. Predictions for the primordial bispectrum (see Section 
14.61 show a similar number of peaks to the power spectrum 
(if we fix two values of £ and vary the third), however, the 
bispectrum fluctuates between positive and ne gative values. 

Therefore, following the method used in ISantos et alJ 

(2003) for the frequentist estimator, we divide the uv plane 
into concentric annuli, of width A. To form a bispectrum 
estimate, three annuli are selected, and we set £1, £2 and £3 
equal to the radii of the midpoints of the annuli. The data 
are binned on a square grid in the uv plane. This binning 
shifts the points very slightly and hence is equivalent to 
relaxing the requirement of exact triangles. We find all the 
possible triangles formed from vectors ui, 1x2 and U3 where 
Ui points to the centre of a cell containing a data point in 
annulus i. We can express our estimator as: 

B(l 1 ,£ 2 ,£ 3 ) = 7" 2 " 3 . (13) 

where F = / 3 |7rcr 2 , the prefactor in equation (II 2t . which 
accounts for the conversion to flux density and the effect of 
the primary beam. The a 2 are used to weight each individual 
triangle of vectors, and are estimated as explained in Section 
14.1.21 Ai denotes annulus i. The visibilities are complex, but 
the resulting bispectrum estimates are real since a(—u) = 
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a* (it), by virtue of the fact that the temperature field is 
purely real. 2 

Our method for estimating the bispectrum from inter- 
ferometer data mir rors that used for the MAXIMA data 
JSantos et alj E003l. However, the MAXIMA data are ob- 
tained in real space and so it is necessary first to perform a 
Fourier transform. We begin with the data in the uv plane, 
so we do not have to concern ourselves with window func- 
tions. However, we have to deal with non-uniform uv plane 
coverage, which results in large variations in the noise on 
each cell. There is also a lot of variation in the number of 
triangles of vectors which form each bispectrum estimate. 



4.1.1 Choice of A 

A natural scale for A, the width of the annuli, is set by the 
width of the aperture function. Increasing A increases the 
number of triangles that form each bispectrum estimate, but 
they correspond to a wider spread of the underlying values of 
£. This loss of At resolution runs the risk of 'washing out' any 
oscillatory signal that may be present. However, if A is too 
small then the estimates have large variances resulting from 
there being few triangles forming each estimate. In addition, 
adjacent estimates would be highly correlated. Therefore we 
chose a value of A which is large enough to encompass the 
width of the aperture function, but with the relative size 
compared to the aperture function slightly smaller for the 
case of the extended array, since it has a larger aperture 
function. Fig. Q shows how particular bispectrum estimates 
vary with the width of the annuli. For the compact array 
we chose A = 16 A, and for the extended array, A = 32.4 
A. This is equivalent to 1.46 and 1.31 times the FWHM of 
the aperture function for the compact and extended arrays 
respectively. 



VSA1E: B(1 000,800,800) 
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Figure 1. Illustration of the variation of a bispectrum estimate 
with A, the width of the annuli, for the extended array field 
VSA1E, with estimated variance calculated according to equa- 
tion 1151 . For small values, there are few triangles and so the 
estimate is large, and fluctuates rapidly. The bispectrum appears 
to converge towards zero as we increase A, and varies smoothly 
with A. We wish to chose the smallest value of A for which the 
bispectrum estimate is changing reasonably slowly, and for which 
adjacent bispectrum estimates are uncorrelatcd. 



(15) 



F 2 E 



In practice, we compute errors and assess the statistical sig- 
nificance of our results using Gaussian simulations which 
properly take account of correlations, however, we can use 
equation 11511 as a way of checking the calculation of the 
weights as given by equation C1J. 



4-1-2 Weighting the bispectrum estimates 

We can achieve a near-optimal estimator by weighting each 
visibility with its inverse variance. In the absence of primary 
beam effects, this estimator would be the optimal cubic esti- 
mator of the bispectrum (see Appendix lA")l . We thus weight 
each individual triangle as ^ 1 where 

^ 1 « 2 »3 = (pi t + [Gt 2 + al 2 ^ [Ct 3 + a 2 u ^j , (14) 

is the product of the variances of the visibilities. Here, Ce i 
is given by ira 2 f 2 C(ui) as shown in equation (1101 . 

We can make a rough estimate of the variance of the bis- 
pectrum by ignoring correlations due to the primary beam, 
and neglecting the fact that some visibilities are used more 
than once in the bispectrum estimate. Under the assumption 
of a Gaussian signal, each triangle is the product of three 
independent visibilities, and nearly all triangles in a given 
estimate B are also independent. It follows that 

2 Although the signal contribution to the visibilities satisfies 
a* (it) = a(—u), the noise is generally uncorrelated between visi- 
bilities. In our analysis we impose the symmetry a*(u) = a(—u) 
on the noise too by averaging the data a(u) and a*(—u) to form 
the visibility at u. 



4.2 Simulations 

4-2.1 Gaussian Simulations 

In order to assess the statistical significance of our results, 
we have performed Gaussian simulations of the VSA data 
from a simulated Gaussian sky. We extract visibilities at 
the same uv positions as our data, and add noise which is 
drawn from a Gaussian distribution with the same variance 
as the noise on the data. For each simulation we compute the 
bispectrum, with the same code as that employed on the real 
data. From the suite of 15000 simulations we estimate the 
variance in each bispectrum estimate, and their covariances 
(which arises from a combination of sample variance and 
instrumental noise). 

The power spectrum for the simulations is drawn from 
a flat ACDM model, w ith the same parameters as model A 
in ISlosar et alJ l|2003|l . using all CMB data. Model A was 
the simplest and had the highest evidence. 

Fig. |2| shows the distribution of some bispectrum es- 
timates obtained from the Gaussian simulations, together 
with the value calculated from the real data. The estimates 
which are formed from more triangles tend to have a smaller 
variance as would be expected. The majority of the distribu- 
tions fit well with a Gaussian curve with the same variance, 
although there are exceptions (which tend to be when some 
of the Is are low) such as the distribution of B(773, 366, 366). 
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Figure 2. Distribution of bispectrum estimates for VSA1F obtained after 15000 simulations, and comparison with 
a Gaussian distribution with the same variance. The measured value from the real data is shown by a vertical line. 



4-2.2 Non-Gaussian Simulations 

It is desirable to test the bispectrum estimator on simulated 
data with a known bispectrum to check for bias. In order to 
do this, we adopt a probability distribution function (pdf) 
derived from the Hilbert space of a linear harmonic oscilla- 
tor, developed bv lRocha et alJ (1200111 . This exercise helps us 
to assess the performance of our code to compute the bis- 
pectrum, and could also be a useful alternative hypothesis 
in testing for non-Gaussianity. We give here a brief account 
of the procedure followed to simulate these non-Gaussian 
CMB maps; a more detailed description will be presented in 
Rocha et al. (in preparation). 

We start by drawing the values of the CMB tempera- 
ture fluctuations, AT(x p ), independently in each real-space 
pixel from our non-Gaussian pdf. Since the pdf is based on 
the wavefunctions of the eigenstates of a linear harmonic 
oscillator, it takes the form of a Gaussian multiplied by 
the square of a (possibly finite) series of Hermite polyno- 
mials where the coefficients a n are used as non-Gaussian 
qualifiers. These amplitudes a n can be writte n as series of 
cumulants dContaldi. Bean~ Magueiidll999h . and can be 
indepe ndently set to zero without mathematical inconsis- 
tency ijRocha et al.ll200ll) . Hence we should regard a n as 
non-perturbative generalisations of cumulants. Let x rep- 
resent a general random variable, within a set of variables 
which are assumed to be independent. The most general 
probability density for the fluctuations in x is thus: 



P(x) = \ip(x)\ 2 = e 



a-nCnHr, 



(16) 



where the H n are the Hermite polynomials, and the quantity 
ctq is the variance associated with the (Gaussian) probability 
distribution for the ground state IV'ol 2 - The C„ are fixed by 
normalising the individual states. The only constraint on the 
amplitudes a n is: 



E 



= i. 



(17) 



This is a simple algebraic expression which can be eliminated 
explicitly by setting ao = yl— • J2T ! a ™l 2 - 

We consider here the situation in which all a n are set 
to zero, except for the real part of Q3 (and consequently 
ao)- The reason for this is that such a quantity reduces 
to the skewness in the perturbative regime. The imaginary 
part of ay is only meaningful in the non-perturbative regime 
(and can be set to zero independently without inconsistency; 
iR.ocha et alJl200lh . Hence we are considering a pdf of the 
form: 




Figure 3. Non-Gaussian pdf given by equation 1181 with ai 
»2 = 0, 03 = 0.2 and ctq = 1. 



P(x) 



-* 2 /(2<^) 



27TO"0 



1 0:3 V 
ao H r=H 3 



s/2, 



o-q 



(18) 



with ao — y/i — of. This allows us to obtain a centred dis- 
tribution with /ii = , where /x n is the nth moment around 
the origin defined as: 



(x n ) 



x n P(x)dx. 



(19) 



The first, second and th ird moments of our pdf are related 
to a-j, and ao as follows dContaldi fc MaeueiidEoOlh : 



in 
Ha 



= 0, 

= cro(l + 6^3) 



= {2<tZ 



3 [<*§(! 



(20) 



Therefore we have generated a centred distribution with a 
fixed variance and skewness. For our purposes here we con- 
sidered the distribution with 03 = 0.2 and 00 = 1 plotted 
in Fig. H 

The space of possible distribution functions is con- 
strained due to restricting the set of parameters to two pa- 
rameters only. This implies that we cannot generate dis- 
tributions with any given variance and skewness. How- 
ever, in general our method can generate higher values of 
the relative skewness (since it can generate any distribu- 
tion) but for that purpose one needs more parameters a n 
jContaldi Ma g ueiidl200l. 

The maps simulated by drawing the pixel values from 
this pdf will be, by construction, statistically isotropic. They 
consist of non-Gaussian white noise with variance given by 
H2- We then Fourier transform the map AT(x p ) to get the 
Fourier modes, a(u), which have variance proportional to 
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m 



Non-Gaussian simulation: B(/ / /) 




Figure 4. A Non-Gaussian simulation of the sky in real space 
obtained using the procedure described in Sec, I4.2."2l 



fi2- We rescale these Fourier coefficients so that the vari- 
ance is given by the correct angular power spectrum Ce- We 
then inverse Fourier transform these coefficients back to real 
space to produce a new signal map, AT(x p ), which now has 
the appropriate covariance matrix. In Fig. 0]we plot one of 
these non-Gaussian maps. This new map can then be used 
as input to simulate the VSA observational strategy and to 
obtain a set of visibilities as observed by the VSA. For our 
purposes we are interested in the simplest case, i.e. with 
no beam convolution and no noise, though these can easily 
be incorporated if we wish. We output the visibilities on a 
square grid, and then compute the bispectrum of these sim- 
ulated visibilities with the same code that we apply to the 
real data. 

By construction, the power spectrum of the non- 
Gaussian map is Ce and the bispectrum is related to the 
skewness and variance of our non-Gaussian pdf by: 



h 1 



—3/2 A>ix V 7 Ce 2 Ce 3 
M 2 



(21) 



where A p i x is the pixel area given by A p i x = L 2 /A r p i x , for a 
small patch of the sky of area L 2 . For a detailed calculation 
see Rocha et al. (in preparation). In Fig. we plot the com- 
puted values of the bispectrum against the predicted ones. 
The agreement shows that our bispectrum calculations are 
indeed correctly obtained. 

We note that from the generalised Baye sian analysis 
using the non- Gaussian likelihood deve loped b y Rocha^^jJ 
( 2001), and applied to the VSA data bv lSavaee et alJ 12004). 
we concluded that VSA data are mostly consistent with zero 
CY3, resulting in no evidence for this type of non-Gaussianity. 



4.3 Testing for non-Gaussianity 

Fig. shows the estimates of B(£,£,£), together with error 
estimates from Gaussian simulations, from the observations 
of the three compact array fields and the three extended ar- 
ray fields in the VSA1 region. Little correlation can be seen 
between the three datasets in each graph. At low values of 
£ noisy estimates are obtained due to the fact that few tri- 
angles can be found at low £. The estimated variance of the 
mean is not dissimilar to the mean itself, indicating that we 
cannot expect to find a significant deviation from zero. The 
bispectrum estimates from the extended array are much less 
noisy than those from the compact array, reflecting the fact 
that the noise on the visibilities is smaller for the extended 
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Figure 5. Bispectrum estimates of the non-Gaussian simulations 
with a\ = a.2 = 0, 03 = 0.2 and ctq = 1 (i.e. with 113 ~ 0.95 for a 
patch of the sky of side 28? 6 and pixel number 7V p j x = 128 2 ), as 
compared to the predicted bispectrum obtained using equation 
1211 , The error on the measured bispectrum are due to the finite 
number of simulations (1865). 



array than for the compact array (the mean noise per binned 
visibility is 4.4-8.8 Jy for the compact array and 0.7-1.1 Jy 
for the extended array). The noise contributes as to the 
bispectrum error. 

Figs [7| and |H| show a sample of the bispectrum values 
computed from the real data, together with the 1 and 2- 
cr values from the Gaussian simulations. The change in the 
variance with multipole is due mainly to the variation in the 
number of triangles of vectors used to form each bispectrum 
estimate. In addition, the noise on each binned visibility 
value varies as a result of the differing numbers of raw vis- 
ibility measurements used to form each binned value. The 
rough estimates of the errors from equation l|15p agree well 
with the l-(j values from the simulations for large values of 
£, tending to be too small by approximately 0-4 per cent. 
However, for small values of £, equation l|15|l generally un- 
derestimates the variance by a greater amount, the most 
extreme example being for 5(366,366,366) where the esti- 
mated variance is only 35 per cent of the value obtained from 
simulations for the case of the VSA1F field. The distribution 
of this particular bispectrum value is very non-Gaussian, as 
shown in Fig. and hence it appears that the correlations 
that we have neglected in equation 11511 are significant in this 
case. Although around 160 different triangles of vectors in 
the uv plane make up the bispectrum estimate, the effective 
number of independent triangles is much smaller. 

For the 17 fields, there are a total of 1839 estimated 
bispectrum values. Of these, 99 were found to be greater 
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Figure 6. Estimated diagonal component of bispectrum from 
the compact (top) and extended (bottom) arrays for the region 
VSA1, with error bars from Gaussian simulations. Some values 
with very large errors at high and low £ have been omitted. 



than 2cr in magnitude, compared with the expected value of 
~84 if the bispectrum estimates were Gaussian distributed 
(which we have found to be the case except for when one of 
the is is small). The VSA3 field had a significant number of 
large bispectrum estimates (18 out of 176 of modulus > 2a, 
6 of which were > 3a), most of the largest of which appear 
to have in common at least one value of £=584. There was 
one bispectrum estimate outside the 4-<r limit, in the field 
VSA2G. In order to proceed further, we need a way of testing 
the data as a whole. In the following section we describe one 
such test that we have applied to the individual fields, as well 
as to sets of bispectrum estimates formed from the weighted 
sum of the estimates for individual fields that are in the 
same region of the sky. 



4-3.1 Kolmogorov-Smirnoff test 

We wish to perform a quantitative comparison between our 
data and the simulations, in order to test the null hypothesis, 
Ho, that the data values were drawn from the same distri- 
b ution as the Gaussian simulations. Following the method 
in lKomatsu et all {2002) , for each bispectrum estimate B a 
we compute the quantity 
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Figure 7. Results from simulations for VSA1 showing variance 
of bispectrum estimates and data values. The top figure shows 
the diagonal component of the bispectrum. 
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Figure 9. Distribution of the bispectrum estimate 
-B(366, 366, 366) for the field VSA1F, and comparison with 
a Gaussian distribution with the same variance. 



P a = 



N ( \B„ 



< \bX sa 



At 



(22) 



where a represents a particular set of {£1,^2,^3} and 
A^T/otai = 15000, the total number of simulations. This gives 
the probability that the magnitude of the bispectrum es- 
timate drawn under Ho is less than the magnitude of the 
bispectrum estimate obtained from the data. If Ho is true, 
then the distribution of P a is uniform on the interval [0,1]. 
If there is a tendency towards a non-zero bispectrum we 
could (naively) expect to obtain more values of P a close to 
1. Fig. llOl shows the resulting cumulative distributions of P a 
obtained for the VSA1 observations, in comparison with the 
straight line expected from a uniform distribution. 
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Figure 8. Bispectrum estimates with variance from simulations for the extended array fields VSA1E and VSA1F. 



VSA1 Compact array 




Fe(P) = 



N(P a < P) 
N 



(24) 



Figure 10. Cumulative distribution of P a for the VSA1 compact 
and extended arrays. 



The Kolmogorov-Smirnoff test operates by finding the 
quantity 

D — max \F e (P) — F(P)\ , (23) 
where F e is the empirical distribution function defined by 



and F(P) is the continuous cumulative distribution function 
for P under Ho- If the P a are independent, identically dis- 
tributed under Ho then the distribution of D is universal, 
depending only on the number of values, N. 

Given D = d, one can then compute the tail probability 
p which is given by 



p(d) = Prob(D >d\H ) 



(25) 



In our case, we must calculate p(d) using the distribution of 
D obtained from simulations, since the assumption that the 
P a are independent is invalid: although the formal covari- 
ance between different bispectrum estimates is zero apart for 
the effects of the primary beam, the same visibility point 
is used to calculate many different bispectrum estimates. 
Fig. 1111 shows the distribution of D obtained from simula- 
tions, compared with the distribution which would be ob- 
tained if the P a were independently drawn from a uniform 
distribution. It can be seen that there is an increase of larger 
values of D. A very similar distribution is obtained from 
noise-only simulations, when there is no correlation between 
any of the visibility points. 

Table^shows the calibrated results. It should be recog- 
nised that this test gives us an indication of how well the 
simulations match the data, which is dependent on a num- 
ber of factors other than simply whether the data are Gaus- 
sian: whether the power spectrum is accurate and the correct 
telescope parameters have been used, and whether the noise 
estimate is correct. For comparison, when we compute the 
values of P a by comparing the VSA1F data with Gaussian 
simulations which have a beam of FWHM 4? 6 rather than 
2°. 05, we obtain p(d) = 0.006. Altering the amplitude of the 
power spectrum by ~10 per cent was found to change only 
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VSA3G: Distribution of D 
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Figure 11. The distribution of D from simulations (solid line) 
as compared with the distribution obtained if the P a were inde- 
pendent (dotted line), together with the measured value from the 
data for the field VSA3G. 



Table 1. d-statistic and corresponding calibrated tail probability 
for the various fields. The number of bispectrum estimates for 
each field is n. The field labelled 'VSA1 Compact' represents the 
weighted average bispectrum from the individual VSA1 compact 
fields, and likewise for the others. 



Field 


n 


d 


p(d) 


VSA1 


171 


0.076 





38 


VSA1A 


106 


0.064 





81 


VSA1B 


106 


0.16 





03 


VSA1 Compact 


172 


0.086 





25 


VSA1E 


75 


0.17 





06 


VSA1F 


75 


0.10 





52 


VSA1G 


75 


0.15 





11 


VSA1 Extended 


75 


0.11 





34 


VSA2 


173 


0.093 





23 


VSA2-OFF 


171 


0.084 





29 


VSA2 Compact 


174 


0.069 





46 


VSA2E 


76 


0.096 





56 


VSA2F 


76 


0.070 





87 


VSA2G 


95 


0.107 





13 


VSA2 Extended 


95 


0.103 





43 


VSA3 


176 


0.118 





05 


VSA3A 


106 


0.068 





76 


VSA3B 


106 


0.050 





97 


VSA3 Compact 


176 


0.084 





23 


VSA3E 


97 


0.107 





45 


VSA3F 


76 


0.090 





62 


VSA3G 


76 


0.176 





05 


VSA3 Extended 


97 


0.111 





39 



slightly the distribution of d, acting to alter p(d) by ~0.01. 
Four out of the 17 individual fields have values of p(d) which 
are fairly low: VSA1B, VSA1E, VSA3 and VSA3G. For the 
total number of fields analysed we would expect an average 
of one to have p(d) < 0.06 by chance. The lowest value of 
p(d) is for the VSA1B data, at 3 per cent. This field was 
observed for a total integration time of only 68 h, in com- 
parison wit h around 200 hours for the majority of compact 
array fields ijTavlor et, al]l2003T ) . and consequently the mean 
error on the visibilities is approximately twice that of the 
other fields. Furthermore, the large value of d is caused by 
an excess of low values of P a , as can be seen in Fig. ED This 
indicates that we should be very hesitant about drawing any 



conclusions about non-Gaussianity in the VSA1B field. The 
large value of d for the VSA3G field is also due to an excess 
of low values of P a , whereas there is an excess of high val- 
ues for the VSA3 and VSA1E fields. The VSA3 data have 
a greater extent in the uv plane than the other two VSA3 
compact array fields (hence the high value of n). On elim- 
inating all the bispectrum values which are not calculated 
for the other two fields, we obtain n — 106, d — 0.149 and 
p(d) = 0.034, so the significance level is increased slightly, 
indicating that the large bispectrum values are not concen- 
trated in the region of large I. The VSA2 compact array 
fields, which have the longest integration times, are consis- 
tent with the simulations. 

By testing the data in this way, we are making no as- 
sumptions about the type of non-Gaussianity that we are 
looking for. This means that our test is very general, but not 
very powerful as it is not tailored for optimal detection of 
any particular non-Gaussian signature. One disadvantage of 
the test is that all of the P a values are given equal consider- 
ation, regardless of the fact that some bispectrum estimates 
are more dominated by noise than others. In the following 
section we consider two tests which are tailored for detecting 
particular types of non-Gaussianity: that arising from point 
sources, and from the time of recombination. 

4.4 Point sources 

44.1 Theory 

A point source of strength s at position xq can be described 
by 

AI(x) = sS(x - x ). (26) 

The contribution this makes to the visibilities is 

Ay(t I ) = sA(5 )e 2 *" '", (27) 

and therefore the contribution to the bispectrum is 

AV(ui)AV(u 2 )AV{u 3 ) = s 3 A, 3 (x )e 2ni£o - iui+U2+U3) 
= s 3 A 3 (x ). (28) 

There will also be cross terms such as V(ui) V(u2)AV(u3) 
but these will have a mean value of zero. In addition, if 
we have several point sources there will be additional cross 
terms introduced. However, on average the contribution to 
the redu ced bispectrum arising from point sources should be 
uniform llKomatsu fc Sr>ergelll200l|) : 

B ps (h,e2,e 3 ) = s ps . (29) 

So we can form an estimate of the point source contribu- 
tion simply by calculating the (weighted) mean value of the 
bispectrum. 

4-4-2 Point sources in the VSA data 

We compared the estimated value of B ps from the extended 
array data before and afte r source subtraction ilTavlor et alJ 
120031 lOrainee et aljEfiolL The results are shown in Table 
[2] and illustrated in Fig. 1121 The mean value of the bispec- 
trum is altered by ~l-6x 10~ 6 /iK 3 when the point sources 
are subtracted. Therefore, although subtracting the sources 
does change the bispectrum estimates, the resulting differ- 
ence is considerably smaller than the standard deviation on 
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Table 2. Mean bispectrum values before and after point source 
subtraction for the three fields in the VSA1 region observed with 
the extended array. 



Field 


BP7(10- 5 /iK 3 ) 
Raw data 


Sources subtracted 


cr/(10-'VK 3 ) 


VSA1E 


1.19 


0.92 


1.19 


VSA1F 
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-1.61 


1.30 
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1.01 
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1.34 
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Figure 12. Diagonal component of the bispectrum for the fields 
VSA1E and VSA1F before and after point source subtraction, 
illustrating the level of change. 



the mean value of the bispectrum which arises from sam- 
ple variance and noise, and so it does not seem possible to 
detect point sources in this way. However, if future VSA ob- 
servations have a much lower noise and probe higher values 
of £ then it may be possible to use the bispectrum to detect 
the presence of point sources. 



4-4-3 Residual point sources 

We can estimate the contribution to the bispectrum from 
residual point sources using t he results of the 15-GHz Ryle 
survey ijWaldram et al.l 12003') . This survey found that the 
differential source count at 15 GHz could be well approxi- 
mated by 



n(S)~K[^-) Jy-V" 1 , 



(30) 



where S is the flux density, f3 = 2.15 and K = 51. Assuming 
a spectral index a — 0.55 we find that, at 34 GHz, K = 30. 

An individual point source contributes to the three- 
point visibility function as described by equation l|28|l . As- 
suming a Poisson distribution of sources, the mean contri- 
bution from unsubtracted sources can be expressed as 



(V(ui)V(u 2 )V(u 3 )) = / A\x) d 2 x J S A n(S) dS 



2 2 St fi 



(31) 



for ui + U2 + us — 0, where S* is the source subtraction 
level. 

For the compact array, taking S, = 80 mjy we ob- 
tain a theoretical contribution from unsubtracted sources of 
~4x 10" 4 Jy 3 (equivalent to 4 x 10"VK 3 ) and for the ex- 
tended array, taking S* = 20 mjy we obtain ~6x 10 _6 Jy 3 
(equivalent to 3 x 10 _7 /xK 3 ). Comparing this with the er- 
ror on the mean value of the bispectrum of ~ 10 _4 /xK 3 for 
the compact array and ~ 10 _5 /iK 3 for the extended array, 
we see that we do not expect the bispectrum to detect the 
presence of the residual sources. 

If we take these values and consider the comparison with 
the predicted bispectrum arising from the coupling between 
the Sun yaev-Zel'dovich effect and w eak lensing effects, as 
given bv lKomatsu fc Spergell i200lT) . we find that the SZ- 
lensing bispectrum will be overwhelmed by residual point 
sources at the source subtraction level of the extended ar- 
ray for £ > 200. At lower values of £ we estimate that the 
SZ-lensing contribution to the bispectrum is approximately 
four orders of magnitude smaller than the error on our bis- 
pectrum estimates. Therefore our data is not sensitive to 
this effect. 



4.5 Simulated point sources 

4-5.1 Noise-only simulations 

As a preliminary test we subtracted a single point source 
from an extended-array simulation with no noise or CMB 
signal. The resulting mean bispectrum 3 was 96 per cent of 
the theoretical value of s 3 A 3 (xo), with all the individual bis- 
pectrum estimates being very similar in value. The discrep- 
ancy is due to the fact that we slightly relax the requirement 
that ui + ii2 + uz = and so the phase in equation l|28|l is 
slightly non-zero. 

On subtracting two point sources (of flux density 1.0 
Jy and 0.89 Jy after attenuation) from the same simula- 
tion, we immediately find that each individual bispectrum 
estimate is very different, varying from -0.048 Jy 3 to -2.4 
Jy 3 , as a result of the cross-terms. The mean bispectrum 
value is -1.3 Jy 3 compared to a theoretical value of -1.7 Jy 3 . 
On subtracting the VSA1E source list, 10 sources in total, 
the mean bispectrum value is -2.95x 10~ 5 Jy 3 in comparison 
with the theoretical value of -3.53x 10~ 5 Jy 3 , and there is less 

3 By 'bispectrum' in this section we mean the value of the three- 
point function as given by equation 1121 since it is natural to work 
in units of Jy when considering point sources. 
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Table 3. The rms noise, mean bispectrum values and standard 
deviations for noise-only simulations with point sources added. 



rms noise / Jy 


<£}/10- 5 Jy 3 


cr/10- 5 Jy 3 




[/10" 5 /iK 3 ] 


[/10"VK 3 ] 


2.2 


8.8 [4.4] 


390 [19] 


0.44 


3.59 [0.18] 


3.1 [0.15] 


0.22 


1.89 [0.093] 


0.39 [0.019] 


0.044 


2.13 [0.11] 


3.1xl0~ 3 [1.5xl0~ 4 ] 


0.0044 


2.22 [0.11] 


3.1xl0~ 6 [1.5X10 -7 ] 



variation between individual bispectrum estimates, suggest- 
ing that with more sources the cross-terms have a greater 
tendency to cancel. 

In order to ascertain the effect of noise on our ability to 
detect point sources we used sets of simulated data with no 
CMB component, just noise, to which we added the point 
sources subtracted from the VSA1E field. Each simulation 
had exactly the same uv and noise template apart from an 
overall scaling factor on the noises. The values of the mean 
bispectra, together with the standard deviation of the mean 
bispectrum obtained from noise-only simulations, are shown 
in Table The mean value of the bispectrum is positive in 
all cases, but a definite detection of point sources is obtained 
only when the rms noise is less than 0.44 Jy. The point 
source bispectrum is completely swamped by noise for the 
case with the highest noise. Fig. 1 1 3l illustrates the bispectrum 
components for each case. 

There were 10 sources in the simulations in total, and 
the sum A 3 (xi)s^ is 3.5 x 10" 5 Jy 3 . This is 1.6 times 

the mean bispectrum value calculated, weighting according 
to the variance from simulations. However, if we use an alter- 
native weighting scheme, according to the number of trian- 
gles of vectors which are combined to form each bispectrum 
estimate, we obtain (for the lowest noise case) a mean bis- 
pectrum value of 2.6 x 10 _5 Jy 3 . The cross-terms introduce a 
lot of variation between the different individual bispectrum 
estimates and hence cause the final values to depend on the 
weighting scheme. 

The rms noise on each visibility point for the VSA1E 
data is ~1 Jy, which is at a level at which the contribution 
from point sources is swamped by noise. In addition, the 
CMB itself will increase further the variance on the mean 
value of the bispectrum. Therefore we cannot expect that 
we will detect the presence of sources in the extended array 
data using this method. 

For comparison, we also performed the Kolmogorov- 
Smirnoff test on the point source simulations. This gave a 
conclusive detection of non-Gaussianity for the cases with a 
rms of 0.0044 Jy and 0.044 Jy but no detection when the 
rms noise was 0.22 Jy [the value of p(d) was 0.42] despite the 
conclusive detection found by computing the mean value of 
the bispectrum. This illustrates how using a tailored statis- 
tic when searching for a particular non-Gaussian signal is 
a much more powerful method of detection than using a 
general statistic. 

4-5.2 CMB simulation with strong sources 

We computed the mean bispectrum from a simulation based 
on the noise and ww-positions of an extended-array field, 
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Figure 13. Bispectrum estimates from noise-only simulations 
with point sources added, with theoretical bispectrum plotted as 
a horizontal line, a represents a particular set of ^2, £3} 



with both CMB and bright sources. The observed sky is 
shown in Fig. 1141 with the point sources clearly visible in 
the map. The mean bispectrum value was 0.3 Jy 3 , a factor 
of 1000 greater than the standard deviation from simula- 
tions with no point sources. The bispectrum is sensitive to 
the point sources as (source strength) 3 and so is useful for 
detecting strong sources but not for weak sources. 



4.6 Primordial Non-Gaussianity 

It has become usual to quantify the primordial non- 
Gaussianity which arises from weak non-linear evolution on 
super-Hubble scales by a single parameter /nl, the non- 
linear coupling parameter: 

TZ(x) = TZg(x) + />n,(Wc(aO - <^gN», (32) 

where lZ(x) is the comoving curvat ure perturbation as de - 
fined according to the convention in lLiddle fc Lvthl fcOOCft . 
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Figure 14. Simulated CMB with point sources, convolved with 
the primary beam. 

and TZg(x) its Ga ussian component 4 llAcauaviva et all2003t 
lMaldacenall200.4l . Even if the fluctuations produced by in- 
flation are perfectly Gaussian, a non-zero /nl will be at 
second-order in perturbatio n theory due to the non-linear 
nature of general relativity (|Bartolo et al J [20041. The cur- 
rent distribution of matter, which is highly non-Gaussian, is 
a result of the non-perturbative non-linear coupling between 
modes, which becomes increasingly significant after recom- 
bination as the perturbations grow. Previous studies have 
not detected any significant primordial non-Gaussianity, but 
have simply put upper limits on it. This is not surprising as 
the predictions from all but the most contrived inflationary 
models are well below the current upper limits. The recent 
results f rom WMAP place an upper limit of 89 on the value 
of /nl llKomatsu et alJ |2Q03i), in comparison with predic- 
tions from slow-roll inflation which give |/nl| ~ lO^-HT 2 , 
apart from the non-Gaussianity introduced by the subse- 
quent evolution of the perturbations. 

We should therefore only expect to place limits on the 
value of /nl using the VSA data, but it is interesting to 
see what limits the data are capable of producing. In ad- 
dition, the techniques developed to estimate /nl can easily 
be applied to any other theoretical bispectrum which has its 
amplitude as the only free parameter. 

We calculated the theoretical bispectrum for 
the case fNL=l using a modif ied version of CAMB 
llLewis. Challinor fc Lasenbvl [20001. The resulting bispec- 
trum has features on the usual acoustic scale of the power 
spectrum, and its diagonal component is shown in Fig. 1151 

In order to relate the theoretical bispectrum to the mea- 
sured three-point visibility function, we performed the inte- 
gral described in Appendix E| which convolves the bispec- 
trum with the aperture function. For each field, we then 
estimate the set of Q a , the theoretical values of the B a for 
/nl = 1 [where a = (^1,^2,^3)] according to equation 

Some authors iKomatsu fc Spergell kopll; ISantos et aLll2003ft 

use a slightly different definition of /nl based on 4>(a;) = ^1Z(x) 
for adiabatic perturbations in radiation domination. This ignores 
neutrino anisotropic stress which produces a five per cent correc- 
tion. 



Theoretical primordial bispectrum for/ NL = 1 
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Figure 15. Diagonal component of the theoretical bispectrum 
for /nl = I- 
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Figure 16. Diagonal component of theoretical measured bispec- 
trum for the VSA2 data with /nl = 1 

using the predicted value of (S(ui)S(u2)S(us)) as calcu- 
lated from the integral. Figure 1161 shows the resulting di- 
agonal component of the theoretical measured bispectrum. 

We estimate the value of /nl from our bispectrum es- 
timates with the estimator 




where the o\ are the variances of the bispectrum estimates 
from Gaussian simulations. We obtain these variances from 
mosaiced simulations to allow us properly to include mo- 
saiced fields. Equation l-'i'-it reduces to the optimal linear 
estimator of /nl from the bispectrum estimates in the limit 
that these are uncorrelated. 

The overall estimate of /nl obtained from the com- 
pact array data is 85, with a standard deviation of 2700. 
For the extended array, we obtain an estimate of -400 with 
a standard deviation of 3500. The 95 per cent confidence 
limits are 5400 and 7000 for the compact and extended ar- 
rays respectively. These contraints are slightly weaker than 
those ob tained from the MA XIMA data using the FFT es- 
timator dSantos et al.ll2003Tl . due probably to the fact that 
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the MAXIMA data extend to lower multipoles with slightly 
higher A£ resolution, although the VSA sky coverage is 
slight ly greater. [The V SA compact array data cover 101 
deg 2 jTavlor et al] |2003f); the MAXIMA data used had an 
area of 60 deg 2 .] This is also why the compact array data 
are better at constraining /nl than the extended array data. 
The bispectrum arising from primordial non-Gaussianity 
falls quickly with I and so extending the measurements to 
higher values of i does not significantly improve the con- 
straints (unless the noise is very low). This is why the 
WMAP data, which are cosmic-variance limited at low £, 
are able to place so much tighter constraints on /nl- 



5 CONCLUSIONS 

Our bispectrum calculations indicate a slightly greater 
discrepancy with the Gaussian simulations for the fields 
VSA1B, VSA1E, VSA3, and VSA3G than we would expect 
if the CMB sky were perfectly Gaussian. However, there is 
little discernible pattern in the way in which the data devi- 
ate from the simulations. A small level of non-Gaussianity 
is to be expected as there will inevitably be a degree of fore- 
ground contamination and unsubtracted point sources, but 
we have shown that in the case of point sources we do not 
anticipate to be able to detect the resulting non-Gaussianity 
with the current level of experimental noise. In the case of 
point sources, the CMB itself acts as a level of noise that 
can only be r educed by taking m easurements at higher I. 
We note that lSavage et all (l2004Tl found some evidence for 
non-Gaussianity in the VSA1 mosaic, which was attributed 
to point sources or contamination from galactic foregrounds. 
It is possible that the excess of large bispectrum values in 
the VSA1E field arises from the same cause. The fact the low 
values of the tail probability p(d) in a Kolmogorov-Smirnoff 
test only appear in isolated individual fields indicate that, if 
the non-Gaussianity which is hinted at is real, it is localised 
in space. 

The limit on the largest scales probed by the VSA 
means that it can only be used to place weak constraints 
on the value of the quadratic non-Gaussianity parameter 
/nl - The data are more suited to detecting non-Gaussianity 
that is present on smaller scales. 
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APPENDIX A: OPTIMAL CUBIC BISPECTRUM ESTIMATORS 

We have estimated the bispectrum using a simple cubic estimator with each visibility weighted by the inverse of its variance 
(see eauation llH^ , This estimator is not optimal, but for interferometer data, where the signal and noise covariances are nearly 
diagonal, it is close to optimal as we shall argue in this appendix. 

We begin by review ing the calculation of the optimal cubic estimator, first given bv lHeavensI (1 19981) and discussed further 
bv ISantos et all l)2003l) . For simplicity, consider real data {a{\ whose three-point function is related to the bispectrum B a 
that we wish to estimate by 

(aidjak) = ^2 QijkB a - (Al) 

Here, a denotes an (ordered) triplet l\ > Li > £3 and Qfj k is totally symmetric in i, j and k. We construct an estimator y a 
which is cubic in the {cii}, 

= H Ei jk a,aja k , (A2) 

ijk 

and which is unbiased, 

(feHX^XX* 5 * =► ^Er jk Qf jk = S a9 . (A3) 

ijk p ijk 

Only the symmeterised part of E°j k enters the estimator y a but we do not impo se total symmetry at this stage as enforcing 
this complicates the variational procedure that follows. Following iHeavensl <ll998t) . we construct the variance of y a under the 
assumption that the data are Gaussian (so the estimator is optimised for only weakly non-Gaussian data), to find 

var(y Q )=^ ^ E^ k E^ j>h > (Cu>Cjj,C kk > + perms), (A4) 

ijk i'j'k' 

where Cy = (cuaj) is the (symmetric) covariance of the data. Unlike IHeavensl dl998fl . we have not imposed symmetry of E^ k , 
so we cannot simplify equation l|A4(l in the manner that he does (see his equation 16). We now vary E"j k to minimise the 
variance of y a , enforcing the constraints on the mean ( equation IA3(l with a set of Lagrange multipliers X a /3, to find 

6^J E(' i j k ^(2C i iiCjjiC kk i + 3CijCk{i'Cji k ')) — A a g(5^. /fc/ = 0. (A5) 

ijk 13 

Note that only the symmetric part E^ k -, enters this expression and so only that part is constrained as expected. We can now 
solve for a symmetric Ef^ k and \ a p to find 

ETjk = Y2 K/3 J2 Q'l fk , (j^C^C-yC^, - 4( - 4 l + N ^ C;,],C y \fiJ k ^j (A6) 
where 

K/3 = ^2 QijkQi'j'k' { C iJ~ C jj> C kki - 4 + N C i' 1 j' C kH C jkj ■ ( A7 ) 

ijk i'j'k' 

In these expressions, N is the number of data poin t s and is typically very large. Our equations 1A6II and lA7i correct the 
results given as equations (21) and (25) in IHeavensl (1998). The error in the latter analysis arises becau se symme t ry of E" h 
is assumed prior to the variation, but is then not enforced during it. However, the terms in error in IHeavensl lll99Sft are 
suppressed by 1/JV for large N, and so the differences from the results derived here are small for N 3> 1. 

For interferometer data, the covariance matrix dj is close to diagonal with correlations limited roughly to the extent of 
the aperture function. Multiplication by C~j reduces to inverse variance weighting the visibilities if the effect of the primary 
beam is neglected. Furthermore, in this limit Q"j k enforces Ui + Uj + u k — 0, so that equation 1A6II . with N ^> 1, reduces to 
the heuristically- weighted estimator of equation l|13|l . In practice, the finite extent of the primary beam makes the heuristic 
estimator somewhat sub-optimal. The necessary developments to compute the optimal estimator are given in Appendix lO 
where we construct Q"j k . However, given the need to perform large- volume Monte-Carlo simulations to assess properly the 
significance of our results, employing the optimal estimator would have been prohibitively slow for the analysis in this paper. 



APPENDIX B: INTERFEROMETER MEASUREMENTS 



The actual values recorded by the VSA are fl ux density measurements in Janskys. The receivers are calibrated primarily 
by observations of Jupiter l|Tavlor et alJl20o3) . with the primary beam A(x) normalised so that A(0) = 1. The intensity 
fluctuations of the CMB can be connected to the temperature fluctuations by 



AI(x,u) 



dB(v,T) 
dT 



AT cmb (x). 



(Bl) 



T=T 



© 2004 RAS, MNRAS 000.IT1IT7I 



16 Sarah Smith et al. 



The visibility seen by the interferometer can be expressed as 

S{u) = j A7(x) A(x) e 27 ""'* d 2 x. (B2) 
The measured visibility V(u) is related to this by 

V(u) = S(u) + N(u), (B3) 
where N(u) is the noise on baseline u. We can relate this to the values of a(u) by 
dB(u,T) 



S(u) = 



T a(u)*A(u) (B4) 

T=T 



dT 

= fa(u)*A(u) (B5) 

where / = 94 x 10 6 Jy sr _1 for v = 34 GHz and a CMB temperature T of 2.726K jMather et al.lll994l) . If we consider the 
variance of the measured signal we obtain 

(V(u)V*(u)) = f 2 Jd 2 ui J d 2 U2 A(u — Ui)A* (u — u 2 ) {a(ui)a* (112)) + cr'i 

= f 2 d 2 ui / d 2 U2 C(u\)8(u\ — U2)A\u — u\)A*{u — U2) + a 2 u 



fC(u) \ d 2 u x \A(u x )\ 2 +al, 



if we make the approximation that C(u) is constant over the width of the aperture function. (We use m to denote |u<|.) Here, 
u 2 , is the variance of the noise, which is statistically independent of the visibility measurement. 
If A(x) = exp(-\x\ 2 /2a 2 ) then we find that 

(V(u)V*(u)> &wa 2 fC(u) + al. (B6) 

Similarly, if we consider 



(V(ui)V(u 2 )V(u 3 )} 



f j d 2 u[ J d 2 u 2 j d 2 u' 3 A(u 1 ~u[)A(u2-U2)A(u s ~u' i )B(u[,u'2,u' 3 )S 2 (u' 1 + u' 2 + u' 3 ) 



f 3 B(u!,u2,u 3 ) / d 2 x / d 2 ui / d 2 U2 / d 2 i*3 



x exp[2-Kix.(ui + u 2 + u' 3 )]A(ui~u' 1 )A(u 2 — U2)A(u 3 — u 3 ) 

f 3 B(u\,U2,u 3 ) / d 2 xA 3 (x) exp[27rii.(ui + 112 + u 3 )], (B7) 



assuming that the bispectrum varies only slowly over the width of the aperture function, and that the noise is Gaussian. 
Taking the above form for A(x) and the case iti + 112 + u 3 — we obtain 

(V(tti)V(ua)V(tt 3 )> w f^wa 2 B(ui,u 2 ,u 3 ). (B8) 
APPENDIX C: VISIBILITY THREE-POINT FUNCTION 

In this appendix we derive an integral expression for the three-point function of the interferometer visibilities. For the special 
case of a Gaussian beam, we are able to reduce the integral further to obtain a closed-form expression for the quantity Q"j k 
introduced in Appendix lAl 

Our starting point is the first equality in equation 1B7II which expresses the three-point function of the visibilities as an 
integral over the bispectrum: 

/oo r 00 p 00 

d 2 ui / d 2 u 2 I d 2 u' 3 B(u[,u2,u 3 )A(u 1 - u' 1 )A(u 2 ~ u' 2 )A(u 3 ~ u' 3 )8 2 (u[ + u' 2 + u' 3 ). (CI) 
-00 J— 00 J — 00 

We can divide the region of integration into six separate regions, u[ > u 2 > u' 3 and its permutations. Since B(u'i, u' 2 , u' 3 ) is 
invariant under permutations of the u' iy as is S 2 (u[ + u' 2 + 1*3), we can permute the Ui instead of the u\ to find 



(S(u 1 )S(u 2 )S(u 3 )} = / d 2 u\d 2 u 2 d 2 u 3 B{u' 1 ,u2,u' 3 ) ^ A(u 1 -u' 1 )A(u2-u' 2 )A(u 3 -u' 3 )S 2 (u' 1 +u' 2 +u' 3 ).(C2) 

J U 1 -^ > > It'? perms 

Transforming to polar coordinates, the delta function restricts the lower bounds of u'2 and by the triangle inequality, hence 



(S{ui)S{u 2 )S(u 3 )} = / uidui / u' 2 du 2 j u 3 dM 3 B(ui,it2,U3) 
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V 



d02 



\ 



t>' 3 ^2 A(ui - u[)A(u 2 — u 2 )A(u 3 - u' 3 )5 2 (u[ + u' 2 + u' 3 ] 

{u 1 ,u 2 ,u 3 } 



.(C3) 



For every u[, u' 2 and u' 3 within the domain of integration, the delta function fixes <j>' 2 = 
ip G [— 7r/2,7r/2] are given by the cosine rule 2uiU 2 cosip = u' 3 2 — u' 2 2 — ii'i 2 , and (f>' 3 = 
Performing the (trivial) integrations over (j> 2 and <f>' 3 , we find that 



i'i + ip, where the two values of 
tfru'-L+u' + n closes the triangle. 



(S(ui)S(u 2 )S(u 3 )) 



Uidu'i / u' 2 du' 2 



u[/2 



u' 3 du 3 



2B(u' 1 ,u' 2 ,u' 3 ) 



K 2 



d0i 



±1/. 



wi)A(u 2 - u' 2 )A(u 3 + u[ + u' 2 ] 



(C4) 



"1 ,1*2 



3} 



with (j> 'o fixed to <j>i + 1/>, an d the additional summation is over the two (±) values of ijj. This result is consistent with equation 
(C2) of lSantos et all <2003h if we transform their result to Fourier space and replace their weight function w(x) by the primary 
beam of the interferometer. 

To make further progress, we specialise to a Gaussian primary beam. As noted in Section^ this is a good approximation 
for the VSA and was assumed for all of the analyses in this paper. For a Gaussian beam with dispersion a, normalised to unity 
at its peak, the aperture function is A(u) = 2na 2 exp(—2ir 2 a 2 v 2 ). If we write u 2 = (u 2 /u'-l)R^,u'-i, where R^, is a right-handed 
rotation through ip, then dot products it, • u' 2 can be written as («2/wi)(RJ Ui) ■ u[. In this manner, the argument of the 
exponential arising from the product of the three aperture functions in equation iCil involves 



(tti - ui) 2 + (u 2 - u 2 ) 2 + (u 3 + ui 



-u 2 ) 2 = E 2 + E' 2 



2u[ ■ U, 



(C5) 



where E 2 = u\ + u 2 + u 3 and E' 2 = u 1 u, 2 i 1*3 , w^v,* ^ _ u, ± u, d vu , 2 / "■i/"-^, 

over <f>i can now be performed to give a modified Bessel function Io(x), and our result for the three-point function reduces to 



u\ 2 + u' 2 2 + u' 3 2 , and the vector U = tti — u 3 + (u' 2 /u' 1 )R, 1 (u 2 — u 3 ). The integration 



{S(ui)S(u 2 )S(u 3 )) = / u'idu'x / u' 2 du 2 

JO Ju[/2 



u 3 du 3 



,2(2n) 4 a 6 B(u' 1 ,u' 2 ,u' 3 ) 



2 W 2 2 



K 2 



■4 2 



4 2 ) 2 



where J7 = |J7|. From this expression, we can easily read off the complex version of Q 

only interested in triangle configurations, u\ + u 2 + u 3 — 0, then it is possible to simplify the product u^U further 



{■Ul,«2.-U3> 

ij k introduced in Appendix^] If we are 



1 TT / 2 / 2 1 2/2. 2/2-, 
U\U =i>\UyU\ +U 2 U 2 + U3M3 J 



-E 2 E' 2 ± AA', 



(ui + u 2 + u 3 = 0), 



(C7) 



where A 



u|)/4 is the area of the triangle formed from u\, u 2 and U3, and A' is the area of the 



triangle formed from the primed vectors. The ± in equation 1U7I arises from the two values of tp. (Whether ± corresponds to 
or =p-0 depends on the orientation of the unprimed triangle, but this is irrelevant for the three-point function since both 
cases are summed over.) 

Finally we note the symmetry properties of the three-point function of the visibilities. For an azimuthally-symmetric 
primary beam, (S(ui)S(u 2 )S(u 3 )) is invariant under reflections, rotations and permutations. For triangle configurations, 
these symmetries ensure that (S(ui)S(u 2 )S(u 3 )} depends only on the lengths ui, u 2 and u 3 , as is apparent in equation 1U7I . 
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